Final Project

library(sf)
library(stringr)
library(dplyr)
library(sp)
library(raster)
library(rgeos)
library(rgdal)
library(rasterVis)
library(ggplot2)
library(terra)
library(sf)
library(fs)
library(here)
library(dplyr)

# TEST IMAGES 
image1B5 <- raster(here("landsat_images","LC08_L2SP_027041_20180829_20200831_02_T1","LC08_L2SP_027041_20180829_20200831_02_T1_SR_B5.TIF"))
image1B4 <- raster(here("landsat_images","LC08_L2SP_027041_20180829_20200831_02_T1","LC08_L2SP_027041_20180829_20200831_02_T1_SR_B4.TIF"))
image1B2 <- raster(here("landsat_images","LC08_L2SP_027041_20180829_20200831_02_T1","LC08_L2SP_027041_20180829_20200831_02_T1_SR_B2.TIF"))
image1B6 <- raster(here("landsat_images","LC08_L2SP_027041_20180829_20200831_02_T1","LC08_L2SP_027041_20180829_20200831_02_T1_SR_B6.TIF"))

image2B5 <- raster(here("landsat_images","LC08_L2SP_027041_20170826_20200903_02_T1","LC08_L2SP_027041_20170826_20200903_02_T1_SR_B5.TIF"))
image2B4 <- raster(here("landsat_images","LC08_L2SP_027041_20170826_20200903_02_T1","LC08_L2SP_027041_20170826_20200903_02_T1_SR_B4.TIF"))
image2B2 <- raster(here("landsat_images","LC08_L2SP_027041_20170826_20200903_02_T1","LC08_L2SP_027041_20170826_20200903_02_T1_SR_B2.TIF"))
image2B6 <- raster(here("landsat_images","LC08_L2SP_027041_20170826_20200903_02_T1","LC08_L2SP_027041_20170826_20200903_02_T1_SR_B6.TIF"))

Pre1B5 <- raster(here("landsat_images","LC08_L2SP_027041_20170911_20200903_02_T1","LC08_L2SP_027041_20170911_20200903_02_T1_SR_B5.TIF"))
Pre1B4 <- raster(here("landsat_images","LC08_L2SP_027041_20170911_20200903_02_T1","LC08_L2SP_027041_20170911_20200903_02_T1_SR_B4.TIF"))
Pre1B2 <- raster(here("landsat_images","LC08_L2SP_027041_20170911_20200903_02_T1","LC08_L2SP_027041_20170911_20200903_02_T1_SR_B2.TIF"))
Pre1B6 <- raster(here("landsat_images","LC08_L2SP_027041_20170911_20200903_02_T1","LC08_L2SP_027041_20170911_20200903_02_T1_SR_B6.TIF"))

#RGB_image <- stack(image1B5,image1B4,image1B2)
#RGB_image <- clamp(RGB_image,0,65535)
#RGB_image <-floor((RGB_image/65535)*255)
#plotRGB(RGB_image,axes=TRUE, main = "False Color: 5,4,2" )

RGB_image2 <- stack(image2B5,image2B4,image2B2)
RGB_image2 <- clamp(RGB_image2,0,65535)
RGB_image2 <-floor((RGB_image2/65535)*255)
plotRGB(RGB_image2,axes=TRUE, main = "False Color: 5,4,2" )

RGB_image <- stack(image1B5,image1B4,image1B2)
RGB_image <- clamp(RGB_image,0,65535)
RGB_image <-floor((RGB_image/65535)*255)
plotRGB(RGB_image,axes=TRUE, main = "False Color Post Hurricane 5,4,2" )

Pre1RGB <- stack(Pre1B5,Pre1B4,Pre1B2)
Pre1RGB <- clamp(Pre1RGB,0,65535)
Pre1RGB <- floor((Pre1RGB/65535)*255)
plotRGB(Pre1RGB,axes=TRUE, main="Fasle Color Pre Hurricane 5,4,2")

#NDVI and NDII analysis

NDVIfun <- function(NIR, red) {
  NDVI <- (NIR - red) / (NIR + red)
  return(NDVI)
}
# Uses SWIR band instead of NIR to look at water absoprtion levels and draw a contrast between them to distinguish rich vegetation
NDIIfun <- function(SWIR,red) {
  NDII <- (SWIR - red) / (SWIR + red)
  return(NDII)
}

#NDII
ndvi <- NDVIfun(Pre1B5,Pre1B4)
ndvi %>%
  plot(.,col=rev(terrain.colors(10)),main ="Landsat NDVI Pre Hurricane Harvey")

NdviPost <- NDVIfun(image1B5,image1B4)
NdviPost %>%
  plot(.,col=rev(terrain.colors(10)),main ="Landsat NDVI Post Hurricane Harvey")

ndvi %>%
  hist(.,breaks =  40, main="NDVI Pre Histogram",xlim =c(-.3,.8))

## Warning in .hist1(x, maxpixels = maxpixels, main = main, plot = plot, ...): 0%
## of the raster cells were used. 100000 values used.

NdviPost%>%
  hist(.,breaks=40,main="NDVI Post Histogram",xlim=c(-.3,.8))

## Warning in .hist1(x, maxpixels = maxpixels, main = main, plot = plot, ...): 0%
## of the raster cells were used. 100000 values used.

#NDII
ndiipre <- NDIIfun(Pre1B6,Pre1B4)
ndiipre %>%
  plot(.,col=rev(terrain.colors(10)),main ="Landsat NDII Pre Hurricane Harvey")

ndiipost <- NDIIfun(image1B6,image1B4)
ndiipost %>%
  plot(.,col=rev(terrain.colors(10)),main ="Landsat NDII Post Hurricane Harvey")

ndiipre %>%
   hist(.,breaks=40,main="NDII Pre Histogram",xlim=c(-.3,.8))

## Warning in .hist1(x, maxpixels = maxpixels, main = main, plot = plot, ...): 0%
## of the raster cells were used. 100000 values used.

ndiipost %>%
  hist(.,breaks=40,main="NDII Post Histogram",xlim=c(-.3,.8))

PreMarB5 <- raster(here("landsat_images","LC08_L2SP_005047_20171003_20200903_02_T1","LC08_L2SP_005047_20171003_20200903_02_T1_SR_B5.TIF"))
PreMarB4 <- raster(here("landsat_images","LC08_L2SP_005047_20171003_20200903_02_T1","LC08_L2SP_005047_20171003_20200903_02_T1_SR_B4.TIF"))
PreMarB2 <- raster(here("landsat_images","LC08_L2SP_005047_20171003_20200903_02_T1","LC08_L2SP_005047_20171003_20200903_02_T1_SR_B2.TIF"))
PreMarB6<-raster(here("landsat_images","LC08_L2SP_005047_20171003_20200903_02_T1","LC08_L2SP_005047_20171003_20200903_02_T1_SR_B6.TIF"))

PostMarB5 <- raster(here("landsat_images","LC08_L2SP_005047_20180904_20200831_02_T1","LC08_L2SP_005047_20180904_20200831_02_T1_SR_B5.TIF"))

PostMarB4 <- raster(here("landsat_images","LC08_L2SP_005047_20180904_20200831_02_T1","LC08_L2SP_005047_20180904_20200831_02_T1_SR_B4.TIF"))

PostMarB2 <- raster(here("landsat_images","LC08_L2SP_005047_20180904_20200831_02_T1","LC08_L2SP_005047_20180904_20200831_02_T1_SR_B2.TIF"))

PostMarB6 <- raster(here("landsat_images","LC08_L2SP_005047_20180904_20200831_02_T1","LC08_L2SP_005047_20180904_20200831_02_T1_SR_B6.TIF"))

Mar1 <- stack(PreMarB5, PreMarB4,PreMarB2)
Mar1 <- clamp(Mar1,0,65535)
Mar1 <-floor((Mar1/65535)*255)
plotRGB(Mar1,axes=TRUE, main = "False Color Pre Maria Hurricane 5,4,2" )

Mar2 <- stack(PostMarB5,PostMarB4,PostMarB2)
Mar2 <- clamp(Mar2,0,65535)
Mar2 <- floor((Mar2/65535)*255)
plotRGB(Mar2,axes=TRUE,main="False Color Post Maria Hurricane 5,4,2")

ndvipreMar <- NDVIfun(PreMarB5,PreMarB4)
ndvipreMar %>%
   plot(.,col=rev(terrain.colors(10)),main ="Landsat NDVI Pre Hurricane Maria")

ndvipostMar <- NDVIfun(PostMarB5,PostMarB4)
ndvipostMar %>%
  plot(.,col=rev(terrain.colors(10)),main="Landsat NDVI Post Hurricane Maria")

ndvipostMar %>%
  hist(.,breaks=40,main="Histrogram NDVI Post Maria")

ndvipreMar%>%
  hist(.,breaks=40,main="Histogram NDVI Pre Maria")

ndiipremar <- NDIIfun(PreMarB6,PreMarB4)
ndiipremar %>%
  plot(.,col=rev(terrain.colors(10)),main="Landsat NDII Pre Hurricane Maria")

ndiipostmar <- NDIIfun(PostMarB6,PostMarB4)
ndiipostmar%>%
  plot(.,col=rev(terrain.colors(10)),main="Landsat NDII Post Hurricane Maria")

ndiipostmar %>%
  hist(.,breaks=40,main="Histogram NDII Pre Maria")

## Warning in .hist1(x, maxpixels = maxpixels, main = main, plot = plot, ...): 0%
## of the raster cells were used. 100000 values used.

ndiipremar%>%
  hist(.,breaks=40,main="Histrogram NDII Post Maria")

PreIan5 <- raster(here("landsat_images","LC08_L2SP_016041_20211011_20211019_02_T1","LC08_L2SP_016041_20211011_20211019_02_T1_SR_B5.TIF"))
PreIan4 <- raster(here("landsat_images","LC08_L2SP_016041_20211011_20211019_02_T1","LC08_L2SP_016041_20211011_20211019_02_T1_SR_B4.TIF"))
PreIan2 <- raster(here("landsat_images","LC08_L2SP_016041_20211011_20211019_02_T1","LC08_L2SP_016041_20211011_20211019_02_T1_SR_B2.TIF"))
PreIan6 <- raster(here("landsat_images","LC08_L2SP_016041_20211011_20211019_02_T1","LC08_L2SP_016041_20211011_20211019_02_T1_SR_B6.TIF"))

PostIan5 <- raster(here("landsat_images","LC09_L2SP_016041_20221107_20230323_02_T1","LC09_L2SP_016041_20221107_20230323_02_T1_SR_B5.TIF"))
PostIan4 <- raster(here("landsat_images","LC09_L2SP_016041_20221107_20230323_02_T1","LC09_L2SP_016041_20221107_20230323_02_T1_SR_B4.TIF"))
PostIan2 <- raster(here("landsat_images","LC09_L2SP_016041_20221107_20230323_02_T1","LC09_L2SP_016041_20221107_20230323_02_T1_SR_B2.TIF"))
PostIan6 <- raster(here("landsat_images","LC09_L2SP_016041_20221107_20230323_02_T1","LC09_L2SP_016041_20221107_20230323_02_T1_SR_B6.TIF"))

Ian1 <- stack(PreIan5, PreIan4,PreIan2)
Ian1 <- clamp(Ian1,0,65535)
Ian1 <-floor((Ian1/65535)*255)
plotRGB(Ian1,axes=TRUE, main = "False Color Pre Ian Hurricane 5,4,2" )

Ian2 <- stack(PostIan5,PostIan4,PostIan2)
Ian2 <- clamp(Ian2,0,65535)
Ian2 <- floor((Ian2/65535)*255)
plotRGB(Ian2,axes=TRUE,main="False Color Post Ian Hurricane 5,4,2")

ndvipreIAN <- NDVIfun(PreIan5,PreIan4)
ndvipreIAN%>%
  plot(.,col=rev(terrain.colors(10)),main="NDVI Pre Hurricane Ian")

ndvipreIAN%>%
  hist(.,breaks=40,main="Histogram NDVI Pre Ian")

ndvipostIAN <- NDVIfun(PostIan5,PostIan4)
ndvipostIAN%>%
  plot(.,col=rev(terrain.colors(10)),main="NDIV Post Hurricane Ian")

ndvipostIAN%>%
  hist(.,breaks=40,main="Histogram NDVI Post IAN")

ndiipreIAN <- NDIIfun(PreIan6,PreIan4)
ndiipreIAN %>%
  plot(.,col=rev(terrain.colors(10)),main="NDII Pre Hurricane Ian")

ndiipreIAN%>%
  hist(.,breaks=40,main="Histogram NDII PRE IAN")

ndiipostIan <- NDIIfun(PostIan6,PostIan4)
ndiipostIan%>%
  plot(.,col=rev(terrain.colors(10)),main="NDII Post Hurricane Ian")

ndiipostIan%>%
  hist(.,breaks=40,main="Histogram NDII POST IAN")

# magian <- ndiipostIan - ndiipreIAN
#magian %>%
#  plot(.,col=rev(terrain.colors(20)))

magian2 <- ndvipostIAN - ndvipreIAN

## Warning in ndvipostIAN - ndvipreIAN: Raster objects have different extents.
## Result for their intersection is returned

#magian2%>%
 # plot(.,col=rev(terrain.colors(20)))



layout_matrix <- matrix(c(1, 2, 3, 4,5,6), nrow = 3, ncol = 2, byrow = TRUE)
layout(layout_matrix)
magmar <- ndvipreMar - ndvipostMar

## Warning in ndvipreMar - ndvipostMar: Raster objects have different extents.
## Result for their intersection is returned

maghar <- NdviPost - ndvi
magianii <- ndiipostIan- ndiipreIAN

## Warning in ndiipostIan - ndiipreIAN: Raster objects have different extents.
## Result for their intersection is returned

magmarii <- ndiipostmar - ndiipremar

## Warning in ndiipostmar - ndiipremar: Raster objects have different extents.
## Result for their intersection is returned

magharii <- ndiipre - ndiipost
hist(magian2, breaks = 40,xlim=c(-0.4,0.6), xlab = "NDVI Difference", ylab = "Number of Pixels", main = "Difference in NDVI")

hist(magmar, breaks = 40,xlim=c(-0.4,0.6), xlab = "NDVI Difference", ylab = "Number of Pixels", main = "Difference in NDVI")

hist(maghar, breaks = 40,xlim=c(-0.4,0.6), xlab = "NDVI Difference", ylab = "Number of Pixels", main = "Difference in NDVI")

hist(magianii, breaks = 40,xlim=c(-0.4,0.6), xlab = "NDII Difference", ylab = "Number of Pixels", main = "Difference in NDII")

hist(magmarii, breaks = 40,xlim=c(-0.4,0.6), xlab = "NDII Difference", ylab = "Number of Pixels", main = "Difference in NDII")

hist(magharii, breaks = 40,xlim=c(-0.4,0.6), xlab = "NDII Difference", ylab = "Number of Pixels", main = "Difference in NDII")